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Ша ИЗТ AL DESCRIPTION OF THE ORBITER 


The space shuttle orbiter is the vanguard of a new 
generation and new concept of svacecraft, having been 
designed to fly as an atmospheric vehicle as well аз to 
operate as a spacecraft in near-earth orbit. In order to 
ШОШО СООП in this dual role, it requires the performance, 


stability, and control systems of both tyres. 


Ворота ilient, attitude control is achieved by 
means of a reaction control system such as those used in 
earlier spacecraft. Whereas earlier vehicles followed an 
unguided re-entry and parachute let-down, the orbiter is to 
be flown to a runway landing and therefore requires the 
conventional aerodynamic roll, pitch, and yaw controls 
ror atmospheric flight. The longitudinal system for 
maintaining the desired pitch angle, pitch rate, and 
At about the pitch axis is the subject of this 


paper. 


The space shuttle orbiter is of conventional design, 
incorporating a low-mounted, highly swept delta wing and 
ENSmatl vertical stabilizer on a thick, rather square body. 
Sea level weight of the vehicle, in re-entry configuration, 


is approximately 186,000 lbs. The wetted area of the wing 





is 2690 square feet, yielding a wing loading of 69 psf. 
Overall length is 122.3 ft, and wing span is 78 ae 
The landing speed of the orbiter, necessarily high because 


of the wing loading, is 190 knots. 


The primary aerodynamic control surfaces are split 
elevons at the trailing edge (to perform the functions of 
ailerons and elevators) and a combination rudder/speed- 
brake at the trailing edge of the vertical stabilizer. 
The rudder splits symmetrically along its vertical hinge 
to produce aerodynamic drag for speed control, hence its 
use as a speed brake. Body flaps are located at the 
trailing edge of the fuselage and, although they were 
originally intended to provide heat shielding for the 
main engines during re-entry, will be used as an active 
longitudinal trim device. Figure 1-1 is a three-view of 
the vehicle and shows the location of the primary control 


Surfaces, 


The elevons are sviit because of wing deflection and 
hinge moment considerations, and are quite large, comprising 
twenty per cent of the total wing area. Their maximum de- 


fkection is 35 degrees up and 20 down; maximum deflection 





1 By comparison, the commercial Douglas DC-9 measures 
119.3 feet in length and has a span of 93.4 feet. 





Figure 1-1 
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Three-view of the Urol ter 





rate is twenty degrees per second. 


The speed brake provides additional longitudinal 
control in some regimes by generating a pitching moment, 


but its primary use i$ for terminal area energy management. 


The elevons are the primary pitch attitude control, 
but since their normal deflection is trailing edge down, 
the body flap is available to reduce the load on them. 
Use of the flap enables the elevons to trail as much as 
possible and minimizes their chances of overheating. 

The body flap is also an aid to longitudinal control at 
high angles of attack where the elevons lose much of 
their effectiveness. Due to rudder blanking, the orbiter 
is in fact directionally unstable in this regime (the 
early re-entry phase). The existence of a lightly damped 
dutch roll mode, however, provides dynamic stability for 


the vehicle, 


The trajectory followed by the vehicle from orbit 
to touchdown is sensitive and complex as with former 
spacecraft and occurs as follows. Follswing de-orbit, 
the reaction control system engines orient the orbiter 


to a normal front-end-forward top-side-up attitude prior 


10 





to reaching re-entry interface at 400,000 feet. The angle 
of attack during re-entry varies between thirty and forty 
degreese When the velocity decreases to 9000 feet per 
second, a slow pitch down maneuver commences which is com- 
ШЕГЕ Бу 70,000 feet; here, the angle of attack has been 
reduced to ten degrees and the velocity to 1500 feet per 
second. Terminal area energy management is then enabled 
to control airspeed and to position the orbiter to arrive 
on the final approach at 12,000 feet, six nautical miles 
from touchdown. Most of the phase prior to the final 
approach position will be flown at a large roll angle and 
at relatively high "g" in order to follow a steeper tra- 
jectory. This results in shortening the range and time to 
touchdown which reduces the total heat load endured by the 
vehicle and extends the lifetime of the thermal protection 


system. 


The flight control system itself has three modes: 
automatic, control stick steering, and direct. Each mode 
may be selected individually for longitudinal and/or 


lateral-directional control. 


The automatic mode consists of hands-off coupling 
of guidance commands through the digital flight control 


Пе tO the control surfaces. Control stick steering 


11 





is a pitch (roll) rate command system for the longitudinal 
(lateral-directional) channel with gain scheduled as a 
function of the dynamic pressure. This mode is further 
compensated for angle of bank and for the pitching moment 
created by the speed brake. The direct system relates 
controller deflections provortionally to the control 
surfaces; this mode will be hopefully unnecessary as the 
vehicle is marginally stable longitudinally, even with a 
forward center of gravity. "The aerodynamic characteristics 
of the orbiter, together with its total devendence on 
MICS for flight control, classify it as a digital 
fly-by-wire, control-configured vehicle."* The entire 
aor oach, from de-orb1t to re-entry to touchdown, will 


normally be flown entirely in the automatic mode. 





2 LCOL en 2 Tuller GON, USAF, “space Shuttle Orbiter Approach 
and Landing Test Program," The Society of Exverimental 
ESL Pilots, 1975 Report, vage 158 
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11. DEFINITIONS AND SYMBOLOGY 


ПЕРІП айт teta cos study, that of automatic control 
has many terms which have a particular meaning when used 
ЕЛ 5 context. So that further references to some of 
these terms is facilitated, the following definitions 
apoly: 

Angle of attack: the angle between the fuselage 
mererence line and the flight path; denoted by ос. 

Block diagram: a graphical representation, such as 
a schematic diagram, of the physical components of a system 
or a set of mathematical equations describing their input/ 
outout relationships. 

Bode plot: a plot of the magnitude of resvonse of a 
system to a Sinusoidal input measured against the fre- 
сасу of that input. 

Characteristic equation: the denominator of a system's 
transfer function set equal to zero. 

Closed-loop: a system in which the input is affected 
or modified by some stage of output. 

Control matrix: describes how the control input affects 
the plant. 

Control system: a set of onysical components related 


in such a manner so as to control or regulate itself or 
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some other system. 

Damping: frictional retardation of the magnitude of 
a deflection (displacement) or rate. 

Damping ratio: the ratio of the amount of damping 
actually present in a system to the amount that would 
critically damp that system. 

Error signal: the actual actuating signal of the system; 
it is generally the difference between the input and feed- 
back signals. 

Feedback: a signal, which is some form of the output 
signal, that is added to or subtracted from the original 
input in order to compare the desired output with that 
actually obtained. 

Eoncing function: see Input. 

Frequency response: a measure of the output of tne 
system as a function of input frequency. 

Gain: the measure of magnitude of an input or feed- 
back signal. 

Initial condition: one of the two types of inputs to 
a system. 

Input: a stimulus applied to a control system in 
order to produce a desired output. 

Linear differential equation: one consisting of the 


sum of linear (first degree) terms; it has the property 
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that the response of a system due to several inputs acting 
at once is equal to the sum of the responses of each input 
acting alone. An equation may be linearized by making 
valid assumptions (approximations) about the system or 

by considering the system to operate in a restricted 
environment where the variations from linearity are 

small (piece-wise linearity). 

Open-loop: a system in which the input is not 
affected or modified by any stage of output. 

Output: the actual response of the system caused 
by the input reacting with the various components of the 
system; lt is not necessarily the response implied by the 
вари. 

Enc: the object or system which is to be controlled. 

Reaction control system: a control system which 
generates rolling, vitching, and/or yawing moments by 
venting jet, rocket, or compressor air exhaust about the 
axis desired. 

Resolvent matrix: denoted by @(s); the matrix which, 
when multiplied by the control vector and the Laplace of 
the input, yields the Laplace of the output of the system. 

Root locus: a method of graphing the roots of a 
Characteristic equation as a function of gain; useful in 


Stability analysis. 


16 


Сару аулам с: the quality of a system that 
determines whether its response to any input will be 
bounded or unbounded. 

Stability, static: the tendency of a system following 
displacement from a rest position, considered positive if 
КОШУ О Лат tendency is to return to the original position. 

State variable: a variable which determines the state 
of the system; the state variables of a dynamic system are 
the smallest set of variables which comvletely determines 
the behavior of the system for any time % 2 0. 

Steady state response: that part of the total resvonse 
which does not approach zero as time approaches infinity. 
Terminal area energy management: a dynamic program 
Mene inputs of orbiter position and velocity to guide 
the vehicle to its final approach vath, minimizing the 

heat load when possible. 

Шиле Constant: inverse of frequency. 

Time-invariant differential equation: a differential 
equation inwhich none of the terms depends explicitly 
upon the independent variable, t. 

Time response: a measure of the output of the system 


function of time. 





3 Only linear, time-invariant differential equations will 
de considered in the development of the equations of 
motion in this study. 
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Transfer function: the ratio of the Laplace transform 
of the output (response) of a system to the Laplace 
transform of the input (forcing function) under the 
Emotion that the initial conditions are Zero. 

iteaaster Mabrix:; the matrix of transfer functions 
for a multiple input/outout system. 

Transient response: that part of the total response 
which approaches zero as time approaches infinity. 

Undamped natural frequency: the frequency at which 
a system would vibrate were it not damped. 

Unit impulse: an input whose magnitude approaches 
infinity and whose duration approaches zero, but whose 


product of magnitude times duration equals one. 


A glossary of symbols follows: 


a rE curve slope, ace 
дес 
b wing span 
© mean aerodynamic chord 
Ст, =, variation in lift with respect to angle of 
a 


attack; generally linear at angles of attack below 


the stall region; coefficient of 
ӘС» 
x ox 
of gravity with respect to angle of attack; must be 


» variation of pitching moment about the center 


negative for positive static stability; coefficient of 
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CMg 


єз] 


Q2 


L 
E 


L,M,N 


pitching moment that remains when lift forces 
equal zero, coefficient of 

EpmeemPxeJYeneyofactor, relative to the elliptical 
ОО and. lirt distribution 

any function, time domain, as f(t) 

any function, Laplace domain, as F(s) 

transfer function, time domain, as g(t) 

transfer function, Laplace domain, as G(s) 

moment of inertia, subscripted as to axis about which 
ie rS taken 

gain 

Laplace transform of, as of [f(t)] 

inverse Laplace transform of, as ШЕСІ 


moments about Х-, Ү-, and 2-ахев, respectively 


ET perturbations in roll, pitch, and yaw rates 


р,а,г 
P,Q,R 


5 


Persurbatıons in roll, pitch, and yaw accelerations 
maces Of roll, pitch, and yaw 


wing area 


u,V,W perturbations in forward, side, and vertical 


velocities 


U,V,w perturbations in forward, side, and vertical 


accelerations 


Oey velocities in forward, side, and vertical directions 


distances along X-, Y-, and Z-axes 
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Х,Ү,2 axis system, right-handed; also force components 
we 7 and Z2 directions 


& angle of attack 


с, 


Dirac delta, unit impulse function d (t) 
deflection, elevon 
damping ratio 


Bien ansle, rate perturbations 


<>. 


pitch angle, relative to earth, positive nose up 


alr denslty 


Oo OP tw 


roll angle, relative to earth, positive right wing 
down (see also resolvent matrix) 


EX frequency 


In order to derive the equatiens of motion of a 
system, a reference frame must be decided upon. The 
most natural of these would be a frame in which the 
axes were aligned with the physical structure of the 
vehicle, one which is called the body-fixed axis system. 
yen a system is depicted in Figure 2-1. Components of 
velocities, forces and moments, and distances relative 
to their respective axial motions are listed in Table 


2-1 9 


In the development of the equations of motion for a 
vehicle, this reference frame has several advantages». 


The rotary inertial properties are constant (assuming 
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applied forces 
motion velocity and moments distance 





forward U X x 
side V I y 
vertical W 7, 7, 
roll P b 
ра Бен 9 M 
yaw R N 
Table 2-1 
X 
U =; 
а 
N 
E 
X-axis 
ONG 
a 
52 М 
Al) „а 
< YA) 
< > 
с? Y-axis 
си 
ПА 
Y 
-axis 
Figure 2-14 


Е Adapted from McRuer, Ashkenas, and Graham, Aircraft 


Dynamics and Automatic Control, Princeton, 1973 





constant mass) so that any derivatives with respect to 
time are zero. The pilot senses forces and moments with 
respect to the body-fixed frame, calculates his reactions 
to vehicular motion and to his changing external environ- 
ment in terms of the body-fixed axes, and generally 
measures his position (attitude), velocities, and accel- 
erations from a body-fixed system. This is especially 
true when he is dealing with primarily relative motion, 


as in alr combat maneuvering. 


An alternative to the body-fixed system is the 
space-fixed system where forces and moments are measured 
relative to an axis system oriented, for example, to a 


flat earth, as in Figure 2-2. 


з--- 
>< 


Figure 2-2 





This system, however, has the very undesirable 
meature of requiring inclusion in the equations of 
motion the variations of moments of inertia of a 
vehicle as the vehicle changes its velocity vector 
and/or attitude. This system would be ideal for an 
aircraft which changed neither, but such is not the 


case in this study. 


Corrections must also be made to the equations of 
motion in the body-fixed axis system, but these corrections 
are limited to inertial forces and accelerations and do 
not involve changing the apparent physical characteristics 
еде vehicle. Further simplifications are enabled by 
this ckoice of systems; one, the small angle assumption 
is made possible leading to the linearization of the 
equations of motion (see Chapter V); two, the symmetry 
of the aircraft across its X-Z plane eliminates several 


terms in the equations. 


The body-fixed axis system is clearly the better 
selection for the development of the equations of motion 


of the vehicle. 


Peder inition of “flight control" is necessary in 
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order to identify those parameters of the vehicle with 
which the study is concerned, The definition is bor- 
rowed from MCRuer, Ashkenas, and Graham : 

"Control: the development and application to a 
vehicle of appropriate forces and moments 
that 

1. establish some equilibrium state of 
vehicle motion (operating point control) 

2. restore a disturbed vehicle to its 
equilibrium (operating point) state 
and/or regulate, within desired limits, 
its departure from operating point 
coude ionese stabilization), " 

Es definition includes the concept of stability and 


effectively excludes the domain of guidanceo 


ШЕНЕП атпстаг: all functions of control (turning, 
climbing, accelerating) are performed by the pilot using 
E Control wheel (or stick), the throttle, and various 
other aerodynamic controls, usually the flap handle, which 
are directly connected to the device which they operate. 
For example, the pilot desires to pull up; he uses some 
degree of back stick, the elevator responds immediately 
and in proportion to the amount of stick travel, and a 


pitch rate is generated. The pilot's environment remains 





? Ibid. 
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relatively constant throughout the flight envelope of the 
airplane; a ten pound pull on the stick will give, for 


example, i.3 g's at a given airspeed in every case. 


Even for higher performance aircraft operating in 
the dynamic environment of air-to-air combat, this is an 
entirely satisfactory system. © Even though a given ele- 
vator deflection will not create the same pitch rate at 
20000 feet as it does at sea level, the pilot quickly 
compensates for this by modifying control stick position. 
However, maintaining a particular pitch angle can better 
be performed by an automatic control system. Using a set 
of accelerometers to sense minute deviations from an 
operating (trim) condition, an automatic control system 
can command the rapid and precise control deflections 
necessary to maintain the vehicle near the desired flight 
Semaition. A pilot attempting to maintain a given pitch 
angle will find his sensory threshold much higher than 
(and his response on the controls much less precise than) 
Meau tomatic system., This type of system is installed on 
many aircraft as an auto-pilot, engaged by the pilot when 


desired, to maintain certain flight conditions such as 





6 Present-day fighters have flight control systems aug- 
mented by 2- or 3-axis stabilization to minimize gust 
response and other oscillations. 
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heading and altitude. 


The development and utilization of control systems 
may be taken one step further. When a vehicle, such as 
the orbiter, is required to fly in vastly different 
environments in the course of a mission, control stick 
inputs to obtain a given maneuver may require very dif- 
ferent control deflections. Variations in dynamic pres- 
sure, Mach number, and altitude create changing values 
of control effectiveness, locations of aerodynamic forces, 
even stability. The automatic control system can use 
inputs such as air density, dynamic pressure, Mach number, 
E existing control deflection to compute precisely how 
much deflection is required to perform a particular maneu- 
ver and how rapidly the deflection should be applied. 

The system may also consider variations in stability under 
different flight conditions. It will then command that 
deflection, measure the response, and make further cor- 
rections if necessary, with far greater accuracy and 


Ac lency than a human pilot. 


This then is the principle by and for which the 
fly-by-wire system was developed, In an environment of 


changing flight conditions and stability where high 
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control precision is required, the automatic system 


significantly outperforms the best of pilots. 
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Ber, MATHEMATICAL MODELS OF LINEAR SYSTEM ELEMENTS 


The Laplace transform is a major tool in the study 
and analysis of automatic control systems. It is the 
procedure by which time domain performance is connected 
to frequency response and results from the correspondence 
between the transfer function and the transient response. 
The transfer function itself is the ratio of the Laplace 
Meeanstorm of à system's output to the transform of its 
ont, assuming that all initial conditions are zero. 

А сгарһіс revresentation of such an input-output rela- 
e lp is provided by the block diagram. In the block 
Eram, of which the transfer function(s) is the heart, 
inputs and feedbacks are represented by their Laplace 
transforms, and the paths which the signals follow are 


depicted. 


Several simplifications in computation occur when 
using the Laplace transform, among which are (1) differ- 
entiation in the time domain is represented by multipli- 
cation by s in the Laplace (frequency) domain; (2) inte- 
eration becomes division by s (both with regard for initial 
mai tions); and (3) convolution, or complex multiplication, 


becomes simple multiplication. 


The equation linking the time domain to the frequency 
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domain 1s 
La со 


lim 
LI] = #6) = 2 |ечва Pe hey at 


ot 


{ 
When the equation of motion for a system can be 
expressed as a linear differential equation, the 
Laplace transform reduces the differential equation 
to an algebraic equation in the Laplace variable, s. 
To illustrate, a damped svring-mass system will be 


used, 


УЛ МЛ лт ае” е” бе” мм 





Figure 3-1 
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Here, my + cy + ky = x(t). Each term is replaced by its 
Laplace transform: 

нд = ғ|ғ70) - 346) - %9 

e | sXG) - y Co) |‏ ج و 

ы = kYG) 

xk) ج‎ X (s) 
When coefficients of like powers of s are collected, the 


equation becomes 


(OY9)]s* - [409 - МӘЛЕ p * euo) - (5) | = X G) 


The appropriate values of the initial conditions are 
entered, and a simple quadratic equation ins has replaced 


a second order differential equation, 


Hecal!ing the definition of a transfer function, all 
initial conditions are set equal to zero, and the ratio 


ШӘ /Х(с) іс found: 
POS Гечев) в + WtG) 7 XO) 
ms a сз + kal v) = KG) 


5-і. = GG) (3-2) 


5 ($) ms" + cs + 





where G(s) is a common symbol for the transfer function. 
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Figure 3-2 is the block diagram of this system where 
ШЕ)! 15 the input or forcing function, Y(s) is the output, 
КОШ С) is the transfer function (also called the plant 
matrix in multiple input/output systems). In this case, 
Meee nas been derived from mass properties of the system, 
the degree of damping, the spring stiffness, and funda- 


mental laws of motion» 


к) с(е) | —=1(8) 


Figure 3-2 


In modeling a system, G(s) will be found in a manner 
analogous to the above. All that remains in order to 
find the response to any input (neglecting feedback for 
ШО теп) іс to multiply the transfer function by the 
Laplace transform of the input. The response may be 
left in the frequency domain, which is the usual case in 
Stability analysis, or transformed back into the time 
domain to study the dynamic characteristics of the 


system, e.g. rise and settling times, maximum overshoot, etc. 


Some of the more common inputs are listed in Table 


3-1 with their time- and frequency-domain representations. 
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The transformations are linear; that is, if F(s) is the 


transform of f(t), then aeF(s) is the transform of asf(t). 


time E 
Leu) X(s) graph 


sinusoidal sin wit 





Table 3-1 


It can be seen in the first column that the unit 
step is the time derivative of the ramp function, and the 
unit impulse (represented by the Dirac delta) is the deri- 
vative of the unit step. If the latter is not obvious, 
consider that the derivative of a function is the slope 
of that function at every point where the function is con- 


tinuous. At t=0, the slove of the unit step approaches 
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infinity, as does the magnitude of the unit impulse. At 
t-0*, the slope of the unit step is zero, as is the magni- 


tude of the impulse function, 


mee all impulse function is an idealization, it 
must be approximated in the real world. When finding the 
Laplace transform of this approximation, any arbitrary 
magnitude may be assigned to the "impulse", as long as 
its duration is small compared to the time constant of 
the system. If the duration of the impulse is taken as 
the inverse of the magnitude, the area under its curve 


will equal one, hence the unit impulse. 


Unless otherwise depicted, successive terms along 
a block diagram path are multiplied together; witness the 
correlation between equation (3-2) and Figure 3-2. This 
is true whether dealing with inputs, successive transfer 
functions, or other operations such as integration, delay, 


adi f ferentiation, Diagramatically, 


X(s) Go (s)| —> Y(s) 
K identical to 
X(s)—»| G4(s)*Ga(s) |-— Y(s) 


E resulting transfer function G(s) = G1(s)°Go(s), or 
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the output Y(s) = X(s)*G,(s)*Go(s). 


A feedback signal, on the other hand, is not mul- 
tiplied but rather employs a summing junction where 
certain signals downstream in the diagram are fed back 
and added to (or subtracted from) the input signal. A 
typical system with one feedback loop is depicted in 
Figure 3-3. The derivation of the complete, or closed 


MOOD) transfer function follows, 


E(s) 





X(s) Y(s) 


Figure 3-3 


E(s) represents an error signal, the difference 
between the input signal and the feedback signal. B(s) 
IHS Called the open loop transfer function; it is the 
input signal after being modified by any transfer func- 
tions and feedback terms up to the summing junction. 

In this case 

B(s) = G(s)eH(s)»*E(s). 
Since Y(s) = G(s)-E(s), 
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and E(s) = X(s) - B(s) » X(s) - G(s)*H(s)*E(s) 
then E(s) + G(s)eH(s)eE(s) = X(s) 
or E(s)|1 A с(в) ен(з) = X(s) 


Ec S= X\s 
(s) 1 + Gls)coH(s 


Substituting the error signal above into Y(s) = G(s)*E(s), 


the transfer function of the system due to feedback be- 
comes 9 _ _Ү\< E G(s 
ше = 7 Í + SETS) ° 


The second order system defined by the damped spring- 
mass system has the unique definition among its properties 


that 


314 


== 252. and = Wn where } = damping ratio 


and Q, = undamped (natural) frequency. If the equation of 
on Is divided through by m, the result is 


Gr gr hy >= Er) 
Оги ü + 3 ہگ‎ 4 + Wn Y = wr x(t) 


where, for convenience, the input function x(t) has ab- 
sorbed a constant 1/k. Taking the Laplace transform as 
before and letting the initial conditions equal zero, 


Y(s) 


Suo. S cu. 


- x 


Ж 22005. 
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If, then, the transfer function of a system may be written 
in the format ъ/(52 + as + b), the physical characteristics 


of the system are immediately evident, 


As will be seen in the analysis of the orbiter, the 
system transfer function may be factored into two such 
second order terms yielding at once the periods and degrees 
of damping of both the phugoid and the short period motions. 
The damping term offers a further clue as to basic or in- 


herent stability of the system. When the roots of the 


quadratics are obtained from 
S = - $o, ы. -1 , the value 


of | determines whether or not those roots lie in the right 


half plane (RHP) according to Table 3-2, 


value of $ roots are 


f< ۵ in RHP 
Дао pure imaginary 
о< 7< ! complex with negative real parts 
= equal, negative, and real 
Y > 4 negative and real 
Table 3-2 


Areeni cance of a root or roots in the right half 


plane will be discussed in the section on stability. 
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ORR STABILITY 


One of the prime indices of a vehicle's perfor- 
mance is the measure of its stability, or how, how much, 
and how quickly it responds to external forces or control 


inputs. 


Two types of stability are generally defined, static 
рис аупатіс. Static stability is the vehicle's tendency 
ENeturn to its original flisht condition after being 
disturbed from that condition. A common representation 
of the three degrees of static stability, positive, 
negative, and neutral, is shown in Figure 4-1. The ball 
possessing positive static stability will return to its 
original position when displaced, the ball with negative 


eerie Stability will continue to move away from its 


(a) (>) (e) 


Figure 4-1 
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original position, and the ball with neutral stability 


will simply remain in its new position. 


The parallel may be drawn with respect to an 
aircraft's longitudinal (pitch) response when displaced 
from trimmed flight? by a control invut or external 
force. Figure 4-2 illustrates the longitudinal forces 
acting on a vehicle where L is the sum of lift forces 
produced by the wing and tail and is located on the 
aircraft so as to create the same total moment about the 
EEuter of mass as would the individual forces. Cm, is 


the pitch moment remaining when the lift forces are zero. 





Figure 4-2 


7 Trimmed flight is defined as that о сонат сов 
where the sum of all moments acting about the lon- 
£ltudinal axis is zero, i.e., Cw = 0. 
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In equilibrium flight, the moment produced by L 
is balanced by Cmo, and no pitch rate is generated. 
Suppose a gust causes a momentary increase in the angle 
lack. If the lift curve slope is vositive (that 
is, a higher angle of attack generates more lift), as 
are all modern conventional airfoils, the increased 
lift will cause the total moment about the center of 
mass to be nose down, or stabilizing. The vehicle will 


Bo return to its original condition. 


If the center of mass should move aft of the lift 
vector, the aircraft could still be stable if Cy, were 
negative (nose down). However, if that gust were to 
hit this aircraft, the increased lift would cause a 
nose UP moment, increasing the angle of attack, further 
increasing the lift, etc. This vehicle then would be 


clearly unstable. 


Three quantities are necessary to this discussion, 
ШИШЕ ТЕ curve slope, the pitching moment coefficient, 
examples of which are shown in Figure 4-3, and Cmg. A 
vehicle must have a positive lift curve slope OE > 


negative vitching moment slope (which means an angle of 


attack above the trim point causes a pitch down and vice 
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En 






Tam point 


ЕСЕ 





AOA 
Figure 4-3 


versa), and a positive Cm, in order to be statically 
stable. Cm, is negative when the aerodynamic center 
of the vehicle lies behind the center of mass, and a 


positive Cm. ls designed into the airframe. 


Given that a vehicle has positive static stabi- 
lity, its dynamic stability may be studied., The 
ball in Figure 4-1(a) is dynamically stable; if dis- 
Maceda from equilibrium, it would return to the bottom, 
probably overshoot, overshoot again coming back (but at 
a smaller magnitude), and eventually come to rest, If 
there were no friction, the system would be dynamically 
neutrally stable; the amount of overshoot would not vary 


from cycle to cycle. If, on the other hand, a magnetic 


8 An exception to this requirement for positive static 
stability is the orbiter itself. Early in the re- 
entry phase, at an angle of attack of 450, the orbiter 
is laterally unstable due to the blanking of the 
vertical stabilizer. However, the presence of a 
Hightly damped Dutch roll gives it positive dynamic 
stability. 
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device were installed under the bowl which was timed to 
attract the ball as it descended and repel it as it as- 
cended, enough velocity could be imparted to the ball 

to finally shoot it over the side. As long as the force 
supplied by the magnet was by any degree greater than 

the forces of friction, the ball would eventually be thrown 
from the bowl; this is an example of a statically stable 


but dynamically unstable system. 


Each degree of stability can be exemplified by an 
rail. A sharp rap on the control stick of a trimmed 
 сгатт causes a pitch oscillation which quickly damps 
out and disappears (positive stability). Most aircraft 
have some degree of Dutch roll, a continuous but mild 
combination of rolling and yawing (neutral stability). 
High performance aircraft, whose stability may be mar- 
ginal under certain conditions, can be susceptible to 
EMO (pilot induced oscillation) where the pilot plays 
the role of the aforementioned magnetic device and 
amplifies the pitching motion of the vehicle, with 
sometimes catastrophic results (negative dynamic sta- 


Баш tv), 


In the more complicated realm of automatic control 
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systems involving feedback, high precision approach and 
landing sequences, etc., more than just "positive" or 
"negative" as regards stability must be known. For one 
reason or another, it may be desireable to utilize a 
nigh gain for a particular feedback signal, but there 
may be a limit as to how high the gain may be before it 
causes the system to become unstable. Certain types of 
inouts, generally frequency related, may cause large 
increases in the magnitude of response (resonance) and 
may require a certain amount of feedback to minimize 


pese excursions, 


Two of the basic tools used for studying the effects 
of feedback on a system are the root locus graph and the 


Bode plot. 


e root locus 15 а locus of roots: the roots are 
those of the characteristic equation of the system. If 
the real varts of these roots are all negative, then the 
tion от Те differential equation describing the 
system will involve negative (decaying) exponents, and 
the system will be stable. This is the rationale be- 
Mme requiring roots to be in the left half plane for 
Stability. Feedback gain can change the location of 


these roots and cause a system to become unstable. 
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An example is in order. Consider the system shown 
in Figure 4-4 where the plant transfer function is 


s-2 
$+1)($+2 


and where the feedback is linear and variable in magni- 


tude. 





Figure 4-4 


The closed loop transfer function becomes 
S-2 
s+1)ís+2) + K(s-2 ; 
and the characteristic equation is 


С БО otic + (222K) = 0. 


The roots, for varying values of K, are given in Table 4-1, 





Fool F1I TOOL #2 


Table 4-1 
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Figure 4-5 is a graph of the root locus of the systeme 
ТЕ may be seen that, as K exceeds the value of 1, one of 


the roots moves into the right half plane. That this 





Figure 4-5 


results in instability is shown by solving the corres- 
ponding differential equation for a value of K>1, Let 


K=2; the characteristic equation becomes 


Converting this equation to the time domain by taking 


the inverse Laplace transform yields 
y+5y-2y=0, 


Bor the solutions to which is y = e0.372t, Since 
y is the output, it may be seen that the output increases 
without bound as the time increases. The system is un- 


stable for this value of feedback gain. 
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The Bode plot depicts how the system will respond 
to a sinusoidal input. Gain and phase margins, which are 
measures of relative stability, are shown directly by the 
plot, as is the magnitude of response at resonance. By 
comparison with templates of the proper scale, the degree 


of damping may be read from the graph. 


The standard representation of the magnitude of 
the response is 20* log |\c( jw)| where jó has replaced s 
in the transfer function expression. This magnitude is 
pot ted against the frequency of the input, which ls on 
a logarithmic scale so that wide frequency ranges may 
be plotted with accuracy. The phase angle, Y , is like- 
wise plotted against the logarithmically scaled frequen- 
су. A stable, time-invariant system, when energized by 
a Sinusoidal input, will have a steady state output of 
the same frequency as the input although the magnitude 


and phase angle will likely be different. 


Since the magnitude of the response is plotted in 
logarithmic form, finding the response curve of a com- 


plicated transfer function such as 


E2 £ 35 + 2 


G mI ISA 
= cJ ToS + ]9s + 13 


a reduced to writing the function as the product of 
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first-order factors, finding the magnitude of each, and 

then adding (not multiplying) the logarithms of the mag- 
nitudes. From this, it is seen that an increase in gain 
merely raises the magnitude curve and does not affect the 


phase angle. 


In the above example, writing the transfer function 


as a function of ju, 


Ge) + Ъцо +2 Qui) + 349 
| 2- —— OT к аге C C = s 
SS o9) « Gu) +13 (0-6) + j au- w’) 


Mewting wo = 0.1 yields 


2.01% 





LAG + уо.» = ç 
A zx — r e a d GC 0.1) = чё 6. 153 o 
Є( |0) = тач + 1.944 да | j | (3.074 

Then 201n|G( j0.1)| = -37.438, the measure of magnitude of 


the response when the input frequency is 0.1 radian per 


second. 


The phase angle is 


11002 En Ma b = 
Re tin aa ¬ Ves a = $.51'- &3$^ « 6.22 . 


Computing the values of the magnitude and phase 
angle over a wide range of frequencies, say from 0.001 


to 10 radians per second, enables the constructién of 
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the Bode diagram. 


ШІ не funcion iS factored into 


= sti )(st+2 
NB) $+3+]2)($+3-]2)($+1 


and each term written as a function of j® , the terms 
become : a A -\ 
jo+1 12-2 [3+3(0-2)] 5+3 6+2) +1. 
The magnitudes at ú = 0.1 аге, DP 

Е 00580.2862 0.273 0.995. 
The logarithms of these magnitudes are 

0.005 0.694 -1.266 -1.298 -0.005. 
which, when summed, total -1.870. Multiplying by 20, as 


before, yields -37.4, the same answer. 


The sum of the phase angles is 
+ 2.669 32,30) = 34,999 - 5.71% = 0,179, 


approximately the same answer as before. 


The factoring method is desired when the Bode 
plot is to be approximated by asymptotes, or straight 
lines which follow the trend of the function. It does 
not consider damping. The former method is entirely 
ШШЕ Степ for plotting the exact function and will be 


the one used in this study. 
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A third type of graph used in the analysis of a 
system is the time response graph, or time history. 
Any forcing function (not just sinusoidal) for which 
a Laplace transform exists can be applied to a system, 
the output calculated, and the output inverted into the 
time domain. This yields a function of magnitude with 
respect to time, the plot of which depicts period(s), 
overshoot(s), rise times, decay times, steady state 
response, and gives an indication of the degree of 
damping in the system. Its drawback is that its degree 
of stability is not easily derived, although it is 


usually apparent whether or not the response is bounded. 


As has been seen, gain will move the locus of 
roots on a root locus plot, sometimes changing the 
stability of a system; gain will move the Bode magnitude 
plot up or down, changing the relative stability. When 
gain is applied to a system, the time history plot 
reflects amplification of the magnitude of the response 
ENNOPoportion to the amount of gain, i.e., twice the gain 
yields twice the magnitude .9 FOF Example, if 5° of 


elevator deflection create 50 lbs of force and 1000 





9 Since we are dealing with linear systems 
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ft-lbs of moment, then 10° of deflection will generate 
Bobs of force and 2000 ft-lbs of moment; the pitch 
rate created by the latter moment will be double that 
created by the former. However, when the gain to be 
varied is in a feedback loop, as described above for 
the root locus and Bode diagrams, the effects on the 
time history are generally unpredictable. The time 
history is therefore of limited value in designing a 
system; its asset is rather in analysis. All three 


types of plots will be used in this study. 
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V. EQUATIONS OF MOTION 


con tinine the study to the longitudinal case, 
only the X and Z forces and the pitching moment, M, 
are of interest. Each of these may have derivatives with 
E D5eet to (1) forward velocity, u, (2) pitch angle, q, 
(3) tical velocity, w, (4) vertical acceleration, w, 
and (5) elevon deflection, бе. These derivatives are 


labelled, respectively, Ks Ха, К Хо,» 2 етс» 


The lower case subscripts refer to perturbations from 


u? 
steady state. 


The linearized equations of motion in state variable 


format, in terms of the derivatives defined above, are: 
= pou ج ا‎ X. = 9 + hou. 
в. = Lau + E = (O, * 2) * іа de 


(5-1) 
(м. А + (Mu ОО 


a 
d 


+ | Ma, + (У, + Ly) Me | E (Me, * То На ) 4. 


ИЕ 
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The state variables are written in the format 
м. 


wh = fF 


and the equations of motion become 


СА X, ри Us E Р 
ve К Ши Да (©, * L.) ° E 
Û 7 м о) (мо Ърм 9 | 
0 2 ° 4 о | [8 
Ks, 
Lee 5, 


ә 


$ (Mg, + Zi Mir) 
LS 
ОИ 1х} = [AM] - (ъ| бе 


where EN is the plant matrix and із) Ponca die neon tool 


тотор «А 


репе justification Is in order for the linearization 
of the system equations of motion. Assumptions are also 
made to simplify the development and/or to define the 
scope of the study. Reference is made to Figure 2-1 for 
sign convention and notation in a body-fixed axis systems 
linear and angular velocities and applied forces and mo- 


ments are identified. 
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The first assumption, which both simplifies the 
development and defines the environment, is that the 
En is fixed in space. In orbital flight, a universal 
inertial frame is necessary due to the spinning of the 
earth on its axis and to the more pronounced effects 
of the gravitational fields of the sun and moon ona 
spacecraft. However, for the relatively short-term 
control analysis of this study, a geo-inertial frame 


OC reference is assumed. 


The airframe is assumed to be rigid. This enables 
the description of vehicular motion as a translation of 
and rotation about the vehicle's center of mass. The 
method of application of the thermal protection layer 
of the orbiter requires an absolute minimum of relative 
motion between and within various sections of the body, 


so this assumption of rigidity is especially valid. 


The mass and mass distribution of the vehicle is 
assumed to be constant. Actually, since no fuel is 
being Ші in producing thrust and no stores expended 
me that part of the trajectory under study, this, 


too, is an especially valid assumption. 


J2 


For body-fixed axes, the inertial equations of motion 
(neglecting gravity) are: 
Te m(U + QW - RV) x 
ЛЕ = m(W+PV-au) $ Z 
ХМ, = AT, + PR Ta) + (PR) La EM 


up 


1 


The gravity force, acting only at the vehicle's center 
of mass, produces no moments but does contribute to the 


Eternal forces; due to gravity, 
AX = ma sm(-8) = - mq sin Ô 
AZ - my cos(-0) соз Ф ma cos Ө сез Ф 


The combination of the inertial and gravitational forces 


м 


on the vehicle yields: 
ра = mm (TS + OW - RV + 45-0) 
Deo с qo or) 


Z 
M 


" 


сие + со тшу) - ROR = ae 


Differentiating these equations yields equations of 


perturbed motion (where, for example, dU = u): 


aX - EE ~ Wee ~ Du Ne - Q. 4 (4058) | 


AL = m li + ap + Po =e au Por 
+ Ca cas Ө, sm Ф, )ф + (m0, cis $.) о | 
AN = à is T (Pe + Rep) (1-12) - (2%. т a. ) “Axe 
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where the quantities subscripted zero are the steady 
state (trimmed) values. To arrive at these equations, 

U was replaced by Ug + u, P by Pg + D, etc, and another 
assumption was made: the disturbances from steady flight 
are small enough so that the sine of an angle may be re- 
placed by the angle itself and the cosine of an angle 
equals one. Any products of perturbations are neglected. 
Me 5та11 perturbation” concept has linearized the 


transcendental equations of the system. 


Because of the symmetry of the airframe in the X-2 
plane, and because there are no aerodynamic asymmetries 
such as propeller slipstream, it is assumed that lon- 
gitudinal forces and moments due to lateral perturba- 
tions about the trim condition are negligible. Pitch 
rate at trim, Qg, is defined to be zero. A further 
Simplification, which now has no effect on the system, 
15 (Пас Ро - Ro = Vo + Фо = zero, This yields 

a = w lš + Wi < 30] 
IL = RS = 9.41 
am = Í 
Since the aerodynamic forces and moments on the 


vehicle are dependent mainly upon the velocity of the 
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vehicle with respect to the air and not upon accelera- 
tion of the air mass, all derivatives with respect to 
acceleration are neglected. The exception is w, an 
important one as it corresponds to the time lag between 
a vertical velocity perturbation at the leading edge 
and its effect at the traiking edge. It would not exist 
ШІ Помпмасвп did not exist, but by this very fact, its 
magnitude cannot be disregarded. The flow is now 


assumed to be quasi-steady. 


Altitude perturbations, being assumed small, allow 


the neglect of air density and temperature perturbations. 


If the body axes are now rotated so that the x-axis 
points into the relative wind, Wo may be set equal to 
zero. This orientation forms a stability axis system. 
(Since small angle assumptions have been made, pertur- 
bations can still be measured in the body~fixed axis 
System.) Applying this rotation of axes to equations 


(5-2) and taking the Laplace transform of each yields 


(5% а = NS +Xo)w < (KS e q os O, do = Ea, Se 
КО + (s- id - 10) w + \CU-2y)s + qm 010 = Zo se 


- Mu - (Mas - М.) сг s ac Мае = Ma Ôe 
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A vector diagram of the forces acting on the ve- 
hicle will aid in the discussion of the derivative terms 
in the above equations; see Figure 5-1. Note that there 


maemo thrust force present. 


L 
x 
ос 3 
© D 
кл 
9 d 

w 

E 

Figure 5-1 


The sums of forces in the x- and z-directions is, 


for small « ы = Dts a = Lax - D 


Var 
7 = - (L + Da) 


EY 


The velocity, V, is assumed to be equal to Ug. 


This assumption is justified by the following treat- 


ment: = 
V Jo + ц)2 + v2 + м2 
y 7o? + 2UU + u£ + vé + we 


1 
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Since u, v, and w are perturbations and are small when 


compared to Ug, their products may be neglected: 
TE US2 4 2Ugu = Uo,/1 + u/Uo = Uo 


Su 
Since АХ = S , then 


24 


30. эр ща ° b 
5. 1 > (тз + sm ne 


= — lA 
За 


where «, is assumed to be zero. Using the definitions 


ү, 2 E and - + p Š, SC, „елеп 
>. AE E + 206.) 
| E vM + c.) 
_ nn a. ab = 
Similarly KZ = u = (EIF t 2р): 
a a.d __1A 
И YA ду = m ou 
Since lift = } = SIC. , then 


453% E SE Me " с.) 


e 
Finally, where the moment = M = > р 9 sz Ln 


А. = n m. + €.) 
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Using the same substitution and differentiation 
techniques as above, the rest of the force derivatives 


Are found. 


The q-derivatives: 








X I POS | Ф Ср Bs 
в” +m 2(%2/25,) 
7 FR PUSE | тб. 
+ ° Кил 2 (5° хл) 
eÓ c Dm 





Му * LCS.) 
The w-derivatives: 


12,5 `e 
шъ O 


M 


0,5 face 
ae == ros >= i Co) 
_ esc | °C. 


Mas 


21. DA 


The w-derivatives: 





Sc ° Co 

Ly * im AL) 
Se обе 

В 

hs * Ци D(2.) 
65 Cm 


Me 7 45 3( 75) 


10The non-dimensionalization procedure is discussed later 
in this chapter. 





The &-derivatives: 
mS 
Ys = - e 


Co 
“о де. 
895 26 
Ode 
= G 


Le = am 
2. 2 
Ms = ed, SE 








Since the forces and moments are proportional to, 
respectively, mass and moment of inertia, the stability 
derivatives are non-dimensionalized by dividing the 
force/moment equations by these quantities; hence, the 


appearance of 1/m and 1/I, in the results. 


Some of the derivatives have particular signifi- 
Васе or are named according to their effect on the 
ШОП of the system. Others may be neglected for 


Zus reasons. A discussion of these items is in 


arder, 


The quantity => is called the speed damping deri- 


Ll Ve, 


The quantity 50 is generally small except at 


transonic speeds, the subject of this paper, 


The quanti ty ES results from aeroelastic effects, 
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which are being neglected, and from compressibility 
ects, which are not. It, too, is largest at 


transonic speeds. 


Cm, is damping in pitch and is negative for a 


q 
stable vehicle. Its primary effect is on short period 


damping. 
е Ха е ° е ° ° 
The quantity =— 15 the static stability deriva- 


tive, the prime indicator of static longitudinal sta- 
bility. Given an increase in angle-of-attack, the 
increased lift on the wing (normally, a tail contri- 
bution is dominant) will cause a negative pitching 
moment (nose down) about the center of mass of the 
vehicle. If this is the case, the center of mass must 
be forward of the aerodynamic center of the wing, and 
the vehicle is statically stable. Considering the 
tallless orbiter, it is expected that Cm, will be small 
compared to a conventional aircraft, becoming somewhat 
more negative (more stable) as the vehicle becomes 
supersonic and the aerodynamic center moves further 


aft on the wing, 


The quantity < is the lift curve slope and is 


positive and generally linear at angles-of-attack 
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below stall, 


Since a change in drag on a horizontal tail is the 
main contributor to a change in Xy, and since this drag 
increment is small anyway compared to the total drag, 


Xw is neglected in further discussion. 


Because of the time lag between a change in angle- 
of-attack at the leading edge and the change in pressure 
кек sc В 
distribution at the trailing edge, edi is produced. 


However, it is à negligible quantity. 


The quantity on will increase short period damping 


when it is negative, as is normally the case. It is 
non-dimensionalized as em 
ACTAS 

The quantities E and = avccsbothesrenificant ın 
tailless vehicles, the first being always positive, the 
second smaller and sometimes negative. Their ratio may 
be likened to L/D. Стр is elevator control effective- 
ness and, as the name implies, is quite importante It 
determines how much pitching moment will be generated 
for a given elevon deflection. It is sensitive to 


center of mass travel, aerodynamic center travel (as 


when a vehicle becomes supersonic), and Сілах" It is 
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a negative quantity for elevons aft of the center of 


gravity. 


Since A = an (Uo) = w/Uo» the terms < and w may 
be interchanged (as may their derivatives) with due 


respect for the constant of proportionality. 


The non-dimensional stability derivatives intro- 
duced on vage 61 will be used wherever practical. 
In accordance with MCRuer's definitions, they are 
called basic derivatives because they do not involve 
any inertial quantities. They are called non-dimensional 
because they involve force and moment coefficients 
rather than the forces and moments themselves. They 
mE useful in that study of the vehicle under different 
flight conditions is most easily correlated. It may be 
seen that, for example, а is the non-dimensional 


form of <, 


mie other type of derivative, of which X, is an 
example, is called a dimensional stability derivative 
parameter. It will be used in the study of airframe 


 апстег functions. 
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VI. COMPUTER METHODS OF ANALYSIS 


The concept of state variable analysis was made 
possible with the advent and refinement of the high 
speed digital computer. The arithmetical operations 
involved in dealing with just a 4x4 matrix are copious; 
they include inverting, factoring, solving the equations 
the matrix represents (which are often differential 
equations), and applying the results to state variable 
procedures for obtaining such things as the resolvent 
matrix, the state transition matrix, and the time res- 
ропве. In addition, most of these results must be com- 
puted for many instances of time, say every 0.1 second 
for 200 seconds, or for a wide range of frequencies in 


order that meaningful data may be plotted and analyzed. 


Before the computer methods can be presented, a 
description of the underlying state variable theory and 


mathematics is necessary. 


Conventional control theory, exemplified by the 
damped spring mass system in chapter III (pp 29 et 


seqo), is adequate for single-input-singl e-output 
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operations but is wholly inadequate for multiple-input- 
multiple-output systems. Each additional parameter 
introduced as an input or output increases the number 


of mathematical operations many-fold. 


The set of state variables must first be specified. 
In accordance with the equations of motion (p 53), the 


ШЕШ М111 ре (u,w,9,9). The state vector is seen to be 


u 
W ; 

e- af» the control vector 1s u = 52%? and the plant 
д 


matrix and control matrix are shown as [A] and {3} res- 


pectively on page 5-2. The state equation for the sys- 


ix] - Lal) + EE . 


The Laplace transform is used to solve this non- 


tem is then 


homogeneous equation of state: 
s Xs) - x(o) = [AIR + ВОС) 
sXG)-WYXG) > KON + BUG) 


ls 1-A] X (5) = xe) + BUG) 
YG) = Us T-AY x (0) m (I-A pu ls) 
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= ще А A 
Since (s I- A) MM. y + 
ЕЕЕ At At 


and n [G1-AY' | = RU - АХ + = + zT о - e 


mien (s) = d£ Ге“ | х (5) + í Ea 85 (s) 


Taking the inverse Laplace, 


t 
t- 
EO. e** x10) ñ En т) 
ò 


where T is a dummy variable., 


buz} dz (6-1) 


wu 


The matrix е^“ is written as (t) and called the 
state transition matrix. Its Laplace counterpart, @(s), 
is called the resolvent matrix. The first term on the 
rieht hand side of equation (6-1) is the time response 
due to initial conditions; the second term is that due 


to the input, Eu. 


Just as in the single-input-single-output system, 
Mere Y(s) = G(s)*X(s), here, X(s) » G(s)*U(s), but G(s) 
is now a matrix and is called the transfer matrix. The 
corresponding closed loop system is depicted in Figure 
6-1. 


Ше) X( s) 








As was developed for the single input/output case, 


В(в) - Н(в)“Х(в) 
ог, dropping the "(s)" for convenience, 
В = НОЕ 
Since X = Go|U - B| 


Go {U - ux] 


chen СНХ = сп 


[I + GoH]X = Gol 


and X =|I + CoH] Cou 


Writing X as CU shows the closed loop transfer matrix 


Ee G(s) = Е + Go(s)*H(s)| +Go(s).11 


In order to solve the equations of motion, a 
computer program will use the adove relationships and 


some form of numerical integration. The programs 


12 


written by Melsa and Jones were used in this study, 


the main two of which were BASMAT and RTRESP. 


'1 Since, in general, [A][3l + Esla], care must be 
exercised in developing the equations by insuring 
that both sides are pre- or post-multiplied by a 
matrix in the same order. 

12 JL. Melsa and S.K. Jones, Computer Programs for 

Computational Assistance in the Study of Linear 

n rol Theory, 2nd edition, MCGraw-Hill, 1973 
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Given a basic matrix, [A], in this case the plant 
matrix, BASMAT produces 
a. the determinant of [А], де+ |А| 
b. the inverse of |А), Lay” 
с. the characteristic polynomial, det(sI-A) 
d, the eigenvalues of EN А: 
Са s Юта оо matrix, (ti) = ae 


«1 
uM ua ssolwentematrix, Ф(з) - (sI-A) 


Once the state transition matrix is known, the 
time resvonse may be found by equation (6-1) whether the 
[put iS in the form of an initial condition(s), a for- 
cing function, or both. The other Melsa-Jones program, 
RTRESP (Rational Time RESPonse), was used for this 


purpose. 


After the mathematical formulation of the time 
response is obtained, 1t is left to graph the results 
to see the time history of the system. In order to 
facilitate this operation, a program named GRAPH was 
written, in BASIC language, for execution on the HP9830 
computer. See Figure 6-2. An explanation of the vari- 


ables and operations in GRAPH follows: 
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MM ER. ROOTS. 


eonj e, 


Bu B LURE RLS 


=й THEM i46 

"out 

BISP "IHPUT RERL ROOT #"1 

igh INPUT PEI] 

134 HEMT 1 | 

140 FOR Т=Е1+: ТО RitR2 STEP Z 

[SG LISP “REAL IMAG PARTS ОЕ СРОК РООТ В 
(ӘН IHMPUT GLIJ«RLIJ 

I BL Ied dec. I3 

[EB RET le BT IJ 

150 HEST Í 

S00 WRITE (152100 

210 FORMAT 450%" 

PRINT “ROOTS OF CHAR EDH ARE” 

PRINT © REAL FART IMAG PART" 
FÜR Ізі TO Fl 

250 WRITE (15, ЗІЙОРГІТ»Й 

266 NEXT 1 

276 PRIHT 


q 


ERU SUE S 
1 Fi е” 
E E Cur fg 


582 ыры 
T 
+ 
bots 


S40 DISP “COEFF OF EEFeUPLDI TY" 
350 ТНРИТ Г 

368 HEXT ! 

STH FOR I=R1+1 TO RisRZ ЗТЕР 2 
Зай R[IJzRBSiRE Il) 

$58 DISP “COEFF DE EXPC"BLIJ"T2sCOSC"RE LI" T3": 
чая [HMPLT CEI] 

418 DISP "COEFF OF EXPC"GLIJ"TOsSIHL" "RET I" T3" 
426 IHPUT CCI*1]j 

438 HEXT I 

418 DISP "MREXIMUM MAGHITUDE OF RESPONSE" 

430 ТНРИТ 59 


Figure 6-2 
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Figure 6-2 (continued) 
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Figure 6-2 (continued) 
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Ci: coefficients of the exponential and sinusoidal 
terms, the sum of which describe the time res- 
ponse 


G9: maximum absolute magnitude of the response; used 
to scale and label the time history plot 


I: counting element 

fas purely real roots 

Qi: real part of complex roots 

Rj: imaginary part of complex roots 
R1: number of purely real roots 

R21 nubes of complex roots 

ШҰ” time; see "X" 

X: horizontal (time) axis 

У: vertical (magnitude) axis 


2: yes(1)/no(0) value; used to answer various 
questions which the program asks of the operator 


es 10-40: dimensions C, P, Q, and R and initializes 
their value at zero 


11 50-310: inputs roots, categorizes them as real or 
complex, and prints results 


11 320-430: inputs coefficients of exponential and 
sinusoidal terms 


11 440-590: scales, draws, and labels the axes of 
the graph 


11 600-650: plots short period magnitude-vs-time 
graph, using subroutine 1080 


11 660-940: plots phugoid graph, using subroutine 1080 
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11 950-1020: allows for labeling 


11 1030-1070: allows drawing of other graphs with 
the same roots 


11 1080-1190: subroutine which computes magnitude 
of response at a given time 


line 1200: END statement 
The values of the coefficients and roots (Ci, Pi, Qi» 


Ri) are obtained from the program PTRESP. 


The text Aircraft Dynamics and Automatic Control, 


referenced previously, contains an input/output analy- 
sis of an early operational Air Force interceptor, the 
F-89 Scorpion. In order to validate the GRAPH program, 
this example was run through the BASMAT and RTRESP ого- 


grams and the results plotted using GRAPH. 


From the given stability derivatives, the plant 
NSEPUX and control vector were constructed in accordance 


with the state equations on page 50 ; for 


altitude = 20000ft 
weight = 305001b 
TAS = 660fps 
Xu = -0, 0097 Zu = -0.0955 Mu = ОТО 

КО = 0.0016 Zw = -1.430 Mw = -0.0235 
Ха, = 0.0 25, = -69.8 Mw = -0.0013 
Ма - -1.920 
Ms. = -26.009 


ща 





the plant matrix [A] was 


20,0097 (00015 0.0 -32.2 
-0.0955 E 660.0 020 
@ 00121 20.0216 -2.778 One 
0.0 0.0 Quo 0.0 


and the control vector was 


-69,8 
-26 з 009 
О 
The input to the system consisted of an elevator 
methection of 0.02 radian commencing at t=0 with a 


duration of one second; initial conditions were equal 


to. The text results are shown in Figure 6-3. 


ШІ огдег to utilize the RTRESP program, a proper 
Laplace transform of the input had to be found. The 
first choice was a multiple of the unit impulse, in 
this case 0.024(t). The Laplace of this input is 
0.02(1/1). However, the Melsa-Jones program requires 
that the denominator of the Laplace of the input be 
one or more orders greater that that of the numerator» 
Since this prerequisite could not be satisfied, it was 
decided to use a step function input and take the deri- 


vative with respect to time of the output, since the 


T? 
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derivative of the step function. The 


system's matrix analysis are given in 
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Figure 6-4. The results of the step function input are 
shown in Figure 6-5, and Table 6-1 shows that the periods 
and damping ratios of the short period and the phugoid 


agree quite closely with the text results. 


Cus ses [ mene Г en 
па fizess | eer [олова | 2.0800” 
f [2923 | oe [Го.от [Гот 


Table 6-1 














The differentiation of the results of the step 


Вас стоп 1п00% у1е1аеа: 


3.718756 


-2.104339 by 


let аі 
0.06274506 


а> = -0,004510971 bo 


u = -0,599148exp(aıt)cos(bıt) 
-0.678176exv(a1 t)sin(b1t) 
ғ0.599209ехр(аә%)сов( 92%) 
%20.06399?7ехр(а2%)віп( 52%) 
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-1,386301 ехр(а1 +) соз(Ъ11+) 
-92,664985exp(aıt)sin(bıt) 

Е leos (lot) 
-0.195852ехр(а>%)віл( 52%) 


= 
th 


-0,520226exp(a1t)cos(b4 t) 
+0,102383exp(aıt)sin(bıt) 
+0,000031ехр(а2%)соз(Ъ2%) 
+0,002468exp(a2t)sin(b>t) 


O 
Il 


+0.039993exp(aıt)cos(bıt) 
-0.117850ехр(а1%)віл( 54%) 
-0,039189exp(a2t)cos(bot) 
-0.,002065exp(a2t)sin(bot) 


D 
II 


Figure 6-6 shows the time history of this system. 
The magnitudes of vertical velocity, pitcn angle, and 
to a lesser degree, pitch rate were decidedly in error. 
meetnermore, the initial condition of zero pitch rate 


Was not reflected in the results. 


The discrepancies are due to the improper descrip- 
tion of the input. Clearly, the effect of a maximum 
elevator deflection acting over a short period of time 
(At« 1) is different from a much smaller deflection 
for one second, hence the magnitude errors. The 
ial condition error is due to the indefinite 
integration of the input. When the input is integrated 
(multiplied by s), the denominator loses a root at 


s=0 and its associated eigenvalue. The time response 
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therefore loses an exponential term, and the initial 


Ponditions are not satisfied. 


The undamped natural frequencies and damping ratios 
are calculated from the roots (eigenvalues) of the 
Characteristic equation by considering the roots in 


mee tor form. 


Мат, 


! 
a а... 
> 
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o 
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Wa = 


м 





where the roots аге ехргеввей ав А -а + jb and were 


calculated by BASMAT. 


The second attempt to duplicate the given input 
function was to evaluate the results of the step function 
(0,02) at t-1 and use these values as initial conditions 
for a Zero-input system commencing at t-1. The values of 
Buesstate variables at t-1, calculated from the RTRESP 


program (see Figure 6-5) were 
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0.9938 
-21.8208 
-0.0035 
-0.0623 


LO 


Ф о жб 


When these values were run through RTRESP as initial 
conditions, with the input function now zero, the 
results were inconsistent with common sense. For 
example, the response of the system a al 
EEuddtion of 0.9938fps of forward velocity perturbation 
and no input yielded a change in U of over 300fps with- 
in ten seconds (the initial velocity of the aircraft was 


but 660fps). 


Using the state transition matrix from BASMAT, a 
program called ICRESP was written to solve the response 


due to an initial condition with the equation 
x(t) = S(t):x(0) . 


The results, shown in Figure 6-7, are much more con- 
Sistent with the earlier figures. In the output por- 
mor of the print-0ut, Al refers to the real part of 
warst root, B? to the imaginary part of the first 
ШОО and likewise for A2 and B2, The matrix E is the 


шШасрах Of coefficients of the state transition matrix, 
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mere the matrix of initial conditions, and X is the 
matrix of coefficients of the output. The time his- 
tory plots of the results are shown in Figure 6-8, 
Since the results were so consistent with the earlier 
program, it appeared as though the reason these figures 
were not coinciding with those of the text was that an 


error existed in one of the two Melsa-Jones programs. 


The results of BASMAT (the determinant, the in- 
verse, the characteristic polynomial, the eigenvalues, 
and the resolvent matrix) were computed by hand and 
mena to be Correct. The program RTRESP calculated the 
same eigenvalues as did BASMAT; its only problem 
appeared to be in the coefficient of the sine term of 
the phugoid root in the forward velocity expression 
when initial conditions are present. The equations 
of state from RTRESP were verified on the HP9830 by 
writing a program called STRESP which calculated the 
response of a system to a step input by using the for- 
mula i 
x(t) = $(+):х(0) + [А] [2($) - nisl.x 
where [A] ist the plant matrix, (BJ the control vector, 


k the gain (magnitude) of the step input, and @(t) the 
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state transition matrix. In this case, x(0) = 0. 


The program and results are shown in Figure 6-9, 


Several other approximations to the desired input 


were tried, the first of which was 


8)) = 1 ер |- т" | hor < < 1 


Transforming into the Laplace domain, 


бо 


sie] = fete е Ге ората 


2 0 
The solution was obtained by the use of Bonnet's 


Theorem (the Second Mean Value theorem) which states 
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mE 


Multiplying by the magnitude of the input, 0.02, 


0.02 $ (+) = —— 


The magnitudes of the results of this input were 


yielded 


too small, the probable reason being that half of the 
Й | 
Вас 51 0п ехр( #4 /т“) lies to the left of the t=0 axis 


CUSI contributing nothing to the invut. 


Ene application or a 0.001 second delay to this 
function would vlace 99% of it to the right of the 
fo axis, countering the problem that existed originally. 
ШОО delay of 0.001 seconds corresponds to multiplying 


-9.0015. 


the Laplace transform by e Approximating 


e-9-001S by the series 


(0.001 NS (0.001) 3n 
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R(s) becomes /f0.000001s? - 0.0005s + 0.0 6.2832 
- 2 
E 52 - 500s + 500,000 


Again, N(s) is of greater order than D(s). 


The correct answer was finally reached when a 
E tion of an earlier idea (using the response at 
Si from a step function as initial conditions) was 
pursued. Іп this case, the two steps were combined into 
one by summing two step functions, one commencing at 
m=O, the other at t-1 and negative in sign. The 
graphical representation and summation is shown in 


Figure 6-14, 





Figure 6-14 


Nine time domain, the result is written 


r(t) = 1(t) - 1(t - тж), where т = 1. 
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The Laplace was easily found and was 


\ "E е: 


606) = 375€ 


Searching for an approximation of e^? which would 
Satisfy the requirement for the order of the numerator 


to be smaller than the order of the denominator led tot? 
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Substituting into R(s) = 0.02(7 - Le?) yielded 


Ris) = 0.04 BE 
S Wao o 55 


The results of the RTRESP computations and the time 
history plots are shown in Figures 6-15 and 6-16. Their 
agreement with the text figures confirmed the accuracy 
of the RTRESP (for other than initial condition inputs) 
and GRAPH programs and provided a valid approximation 


of a finite impulse-type input. 
120zata, Modern Control Engineering, Prentice-Hall, 
1970, p Да 
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VII. LONGITUDINAL FEEDBACK 


The primary purpose in using feedback in a control 
system is to reduce the system's sensitivity to external 
perturbations. In the open loop system, the area of 
main concern is calibration; in the closed loop sys- 
tem, it is stability. Use of feedback can enhance the 
rapidity of response, reduce errors of response, and in- 
crease damping to satisfactory levels; it can also over- 
correct errors or decrease damping to the point where an 
otherwise stable system becomes unstable. It is of in- 
terest, then, to be able to "measure" the stability of 


a system. 


The stability of a system may be ascertained by 
an analysis of the characteristic equation of the system, 
ШЕТ 15, the denominator of its transfer function. 
Famply stated, if the real part of every root of the 
equation is negative in sign, the system is stable. 
That this is so is manifested by the EE between 
the roots of the characteristic equation and the solution 


of the differential equation which it represents. 
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For a multiple input/output system, the feedback 
effects on the transfer functions must be dealt with by 


Matrix manipulation 
The solution to the Бен of differential equations 
x = Ux pu 
was seen to be, in the Laplace domain, 


x(s) = Bra] ho - [$ (=) ъјасә). 


E oven loop transfer function was then 


EONS‏ ا ج 


u(s) 


where @(s) is the resolvent matrix. Ina 4x4 system, 


MN ransfer function is written 
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Ф, Фа Ф ly p. РА Qij b; / 5 





Lv, Qu Ф, © p. 2 Paj b; / A 
(О (5 )= ©», Loa Фу ©, b; I y OR b;/A 
JW 
Qu Pur © Pe 
чъ + bu > Фа; b; /A 
А 
J: k. Ë 


The last quantity will be written 


5 
Сі, where A = det [sI-A] 


Applying the feedback signal to a particular parameter, 


for example, Хз, yields 
Н = lo 0 n ol 
where the block diagram is 


x (S) 





(The solid arrows are used to signify a multiple input/ 


Output system.) The closed loop transfer function is 
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Cals) 
G(s) Е 1 + GoH 


In matrix format, ES) = Е + во | |со(е)) . 


To find this matrix product, the matrices were taxen 


term by term: 





1 0 hG4 0 
i O 1 nG2 0 
|І + оң) -|0 0 163 0 
0 O ПО ] 
e 0 El 0 
0 1+nG -hG 0 
[г + со = | о оз: о 
0 0 “hey, 1+hG3 
А 
Where A = det (I + GoH| = 1 + hG3. 
Finally, 
6, L b; V Z b; Ф\; 
G. Z * Pa; 2 xg 
©. b: Y a 
^ open ge ЖШ, by Pes 
(4e) (eub: — = ——— = алайы ы 
А 
Z b; oo 
2 у; Paj 
2 bj Фа) 
E 2 офу, 
A+ WN 2 bj 9; 
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The closed loop transfer function is then 


0 + 0 Ф. + as + bs Gy 
5,4, 1 LN ar + b3 Ф.з + by Pig 
b, 3, * 5,94, + b, O33 EU by Pry 


e uon her + by Ay + bu Pay 





teo] - 


A + NON x б, фу, + Vo + bu Pas) 


In the general form, the first subscript of Y ín the 
denominator will vary as the feedback parameter is 
varied. Regardless of the form or magnitude of h, only 
the denominator is affected; the mn 1S 4 function 
solely of the openloop resolvent matrix, ®(s), and the 


КОПГО! vector, {b}. 


The most important characteristic of the transfer 
function is that the frequency response may be obtained 
by letting s - jO , where o is the applied sinusoidal 
forcing function. That this is so is shown by assuming 
E the output contains a sinusoidal component of the 


same frequency as the input in the form 
АЕ (4) ви | + + Q (lv j| 


where A is a constant, R(o) is the magnitude, and oe) 
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is the phase angle between the input and the response. 
For an input of AsinWt, the Laplace of which is 


AQ 
сш" 


mae output is 


E), CS | 
у(=) = С(5)г(в) < m EE 4 J; =) 


Using partial fraction expansion: 


ym Де 
y(s) = Т 2 s CIS MESAS 


Assuming the transient solution ls stable allows the 


transient terms to be dropped, so 











Au Els) - Ае(-4ч) 
| $590 $> -ju i 2) 
) 
E I AG (ju) 
г %) A) 


Solving for y(t) by taking the inverse Laplace transform: 


-Ae(-j9) ut ; at 
y(t) = л» + А6 0%) е 


-1 
x x e (7-1) 
Since G(j™) is complex, write it as 


(30) = x(w) + jYlw) and G(-jo) = x(o) - jY(9) 


6. 
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where |) j») | = |X2 (5) + У2(%) and argG( je) = tan- s : 


Substitute into equation (7-1): 


Am = ХЫ je) e = [Xa >] е^“ 


= 
- 


AX GJ D e) = en Ба : 4” 
2) 2) 


A X(o) $^ dt + AYG o9 wt 


м 


AL Xx) + “Y Да [ut + a(w) | 


AMO 
where x = I% X (o) 


This is the same form as the assumed input, so that 
R(s) becomes | Х2(ы) + Y2(w) = | c( jw)| and the phase 


angle y (3) becomes tan”!Y(w)/X(w) = argc(jv). 


It becomes a simple matter to apply a feedback to 
a plant пав корат the transier functions in each of 
the variables, and plot the magnitude of the responses 


to a Sinusoidal input against the frequency of that 
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ul. This plot is termed a Bode ( B6' -dee) plot, named 
for its developer, H.W. Bode. One of the important 
Mees of the plot is in measuring the stability of a 


System, 


In order to facilitate construction of a Bode plot, 
and in order to avoid relying uvon asymptotic approxi- 
mations of magmitude and phase angle, a program named 
BODE was written for the HP9830 computer/plotter. Using, 
as input, the stability derivatives and the type and mag- 
ude of feedback, the program prints out the resulting 
transfer functions (in the Laplace domain), repeats the 
feedback inputs, and, for given limits of frequency and 
magnitude, constructs the Bode diagram on semi-log paper. 
The feedback form may be linear, integration (1/s), ог. 


differentiation (s). 
The program, listed in Figure 7-2, is explained: 


line 10 dimensions open and closed loop matrices 


11 20-60 defines variables in terms of stability 


derivatives; 
XU = Xu 
XW = Xy 
XDE = Xk 


MWDOT = W , etc. 


line 70 iterates program for Xj, X2, хз, and xi 
(u,w,q0,8) outputs 


108 





AS 
p | Ja ° |, |: e “TE j: 
x TINTE 
> FIR Mts MME 


Rete te Б ee 
Е ЕРРЕТИ: ИА 
= i 


a a 


una E a 
NI E eet = 


: 5 
|°. . 


- - 
3 
- .. am 
— 


... 


дръж за е: Е 


-d 
= 
— --- wu 
E = үз 
e E = 


i 
' 
E 


m л 9 ы Е) 
liklz j bE j Fe 


. ere 
+ "emak" 


Ш. 

Alak, Eu 
E ud, На i I, 
i] 


И Eaa 
ae ee a EE 


ve tuners 


a: 
-i 


.. REL لي اسم‎ ie ы — s $ 


22 
Ж 
. . 
.. an: 
2 о 
ач P cd omm | 
DS E 
CO bed Cp 3% = 
ДЕН 5 mi cab iu 
4د‎ К - 1. , А 
| А e Me | E г Ра |1 ы! de 
КОО 4: 
ша 1. Ро 1 ы! 
. .. -— - “. 
Ei - gs “e e ےا‎ 
1. a u 1 
NS. 


E bec bb db ` 
г uu uni Ë 4.” 
г грее РН 
КО 15-1 22 
I E =. 
E. ei 1 г: Ке ч 1 


г” 
A 


E E ИЕ 


--- 


Е (id 1 i Dm 


i 
yore 
1 


i me гт 


t 
= 
т 
un 


el de e 


Кой HEK] 
Кш ПГ]: 
MA sd dere el seh 
сас HE 44 2 шаа 
Ail HE qa 1 Jess 
ED 4]1e-Ma-Zz-4l 
за HL ad done lea MAPA ASA MA AE. 
ив Ala 2 dealer inner 
ala HE. | 
ч ' 


ep ict 
Cuin quor 
EE 


rT 
= 
-+ 


, 
М EXÎ 


] 
SEeMS-MES 21+ NLS CNS +22 9-M4-M Les 

I E E МЕКО НОЕ 
“i 

] 


оры 
Иык ы 


A к хеее А вид РИ РР РР За 


A A E 


ІМ» me 


i m ish i: f HL 
м 
Ера ее ра 


НА 
ا سے‎ 
ұз 

— 


г 
-- 


Е! 


1 
UI 
[ei ee PORTO b 


АНЕ ба ей 
| 


„кж“, 5 
... . » 
. 4 4.7 ' 


жі аға 
ы 


z-H&s4l-MesZIeMISHIAMÉe23 


A a a E A Нох #29 


AS 


м“ 
un 
zl Pet 


E 
fet tat iat 


` 


Figure 7-2 
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ШІ 80-110 format 


11 120-270 assigns values to variables, initializes 
T3 (phase angle, used for comparison) 


11 280-330 initializes numerator and denominator 
coefficients to zero 


11 340-540 computes values of numerator and denomi- 
nator coefficients 


11 550-1070 inputs feedback type and magnitude, com- 
putes and prints closed loop numerator and 
denominator coefficients 

11 1080-1640 inputs frequency and magnitude limits 

(of response), labels graph, and plots 
magnitude and phase angle against applied 
frequency 


11 1650-1900 subroutine for scaling and labelling 
magnitude gravh 


11 1910-2070 subroutine for scaling and labelling 
phase angle grapn 


11 2080-2280 subroutine which computes phase angle 
correct for sicn and quadrant 


line 2290 END statement 


To verify the program, the MCRuer example was again 
used. The text results are shown in Figure 7-3. The 
BODE results are shown in Figures 7-4, 7-5, and 7-6 for 
the three parameters common to both programs, д, х, апа 


uo 


The comparison of the « values of the text and the 


graph of w from BODE needs some explanation. The text's 
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magnitude for « at low frequency is “4,0; the results 
of BODE showed the magnitude of w at low frequency to be 


“61, a 57 decibel difference. For U, = 660, 


> ا‎ > OOO 
log(w) = log(660) + log (a). 


Since log(660) = 56.4, the 57 db difference is correct. 


шаше 7=2 compares the results of the BODE program's 


meester function coefficients with those of the text. 


NUMERATORS 
css = 791, 51085 = 111%.4512 


0.110352 - 7998 - 1151 


69.853 + 26675.783з2 + 258.7185з + 79.9357 


69.852 + 1735052 + 168,65 + 80.2 





BODE 26.00926s* + 35.93499s + 0.35011 


text 26.1252 + 36.15 + 0.351 


DENOMINATOR 
BODE s* + 3.359783 + 26.6754s2 + 0.26271s + 0.07198 
text st + 4.2253 + 18.2852 + 0.20975 + 0.0724 


Table 7-2 


I 





Inasmuch as the text figures appeared to be slide-rule 
accuracy, the results of BODE were considered more 
accurate. The BODE graphs, which were computed with 
zero feedback, compared very favorably with the text, 


and the program may be considered verified. 


The transfer functions may also be written as 
ШИС оп of the stability derivatives, and, in this, 
way, the effects of various types of feedback may be 
predicted. Writing the equations of motion (equations 


5-1) in the Laplace domain: 


(s-X,, Ju -XwW -Xqq + gð = X 
-Zyu + (s-Zy)w -Uga = 24, 
-M U -Мим + (з-Ма )а = Mg, 

-а + 50 = 0 


ШЕ ОТ Сгатег!с rule allows for the solution of each 
variable in terms of the stability derivatives, That 


is: 


Хе -Xur “Ky % 

1. 5-1: ла о 

Қы "Mur $-М 0 

ul) _ Ó ° “| S 
SS А 
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S- Xu Aus - К. 4 

-7, s- L... -Us 0 

where A = Mu Nu Мао 
[о © zi 5 


The transfer functions for w(s)/de(s), a(s)/de(s), and 
D(s)/óe(s) are similarly found. Compiling the coef- 


ficients for each of the numerators yields Table 7-3. 


The time history response may be obtained by sub- 
Titling the error signal for the input term to the 
Ent matrix. In other words, where, for the open loop 


system, 


k(t) = [alx + (и, 


ME ihe closed loop system 


ae) 





ШЕТ - кх, and x(t) = [А]х + {o} (r - Ux). 


Rewriting, . 
x(t) = (A - dk%)x + br. 
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This is the new vlant matrix from which the eigenvalues 
may be calculated and roots plotted. As ky, ls varied, the 
roots move,.and this movement is the critical factor in 

Ш пе, gaining, or improving stability. Plotting 

the roots also shows how variations in k affect the 
undamped natural frequencies and damping ratios and gives 


a rough idea of the sensitivity of the system. 


In order that the roots of a given system subject 
to feedback could be quickly ascertained, a program 
named FDBK was written for the HP9830. Using as inputs 
the resolvent matrix of the open loop system, the control 
MEC tor ío], the coefficients of the open loop character- 
istic equation, and the type and magnitude (gain) of the 
feedback, the program computes La - vx], finds and orints 
the new roots, and calculates the resulting frequencies, 
periods, and damping ratios. See Figure 7-7 for a listing 
of FDBK. The REMark statements in the listing provide an 


explanation of the procedures used. 


Figure 7-8 is an example of the program's output. 
In this case, the M“Ruer example was again used with no 
feedback applied; the results were in agreement with the 


example. 
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Figure 7-9 (a,b,c,d) shows how the roots vary 
as the feedback gain is varied and as different para- 
meters are fed back. Complex conjugates are not shown 
where their omission would be confusing. These results 


E tabulated for more careful scrutiny in Table 7-4, 


Figure 7-10(a) depicts the time history of the 
MCRuer F-89 aircraft subject to an initial condition 
memvertical velocity = 10 fos. No control inputs are 
postulated. The parameter chosen to disvlay system 
response is pitch angle; it was chosen because of its 


good definition of short period and phugoid motion. 


Feedback of horizontal velocity error can sub- 
stantially increase phugoid damping before the short 
period mode is much affected. However, a further 
increase in the gain of this error signal will drive 
the short period to instability. The sensitivity of 
the system in this regime is shown by Figure 7-10(b); 
Mere, With k, = 0.02, the short period motion is less 
than critically damped and oscillates several times 
before dying out. The damping effect on the phugoid 


is readily avparent. 


Feedback of vertical velocity (or angle of attack), 
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Feedback parameter: HORIZONTAL VELOCITY (u) 


cain 


Feedback parameter: 


gain 
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Feedback parameter: 
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if the gain is small, will increase the frequency and 
damping ratio of the short period mode. Since the angle 
attack Signal is 1/660th of the vertical velocity 
Signal, it would be a better parameter from which to 
obtain an error signal. Increasing the gain of this 
feedback will cause a third mode to go unstable; however, 
its period is so great that no problem would be antici- 
Bed in controlling it manually. A further increase 

in gain will cause a fourth mode of significantly shorter 
period to go unstable. Since its period decreases rapidly 
with an increase in gain, it may be assumed that values 
of gain in this region are to be avoided. Figure 7-10 
(c) shows a time history of pitch angle versus time with 


ky = 0.0001 апа мо = 10fps. 


HeeGback Of pitch rate has almost no effect upon 
e рһисоіа motion but appreciably increases the short 
period damping. This very desirable feature makes this 
parameter useful in speed, altitude, and attitude hold 
modes for auto-pilots. Figure 7-10(d) shows a time 
history of pitch angle against time with Ка Ol and 


Wo = 10fpS. 
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Pitch angle feedback is also an excellent method 
of improving longitudinal response. In this case, the 
phugoid damping is significantly increased at the expense 
of short period damping. If short period damping is mar- 
ginal to begin with, use of this feedback may cause an 
instability in that mode. Figure 7-10(e) shows the 
SIS Of Kg = 0.05 on the time history of the pitch 


angle perturbations with wo = 10fps. 


Use of both pitch rate and pitch angle feedback 
together should provide an ideal system with maximum 
increase in damping ratios and minimum increase in 
frequency. USing the same feedback gains as above, the 


response is shown in Figure 7-10(f). 


The equations of state for each of the above 
examples were obtained by entering BASMAT with the 
Mecitied plant matrix [А - ъкт | and multiplying the 
Meet cine state transition matrix by the initial condition 


vector in accordance with the equation 


x(t) = &(t)x(0) where x(0) = 


One further conclusion was reached in this treatment 
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of feedback: the two parameters which are most affected 
in phugoid motion, forward velocity and pitch angle, were 
the two whose feedback most affected that motion. Their 
effect on the short period was minimal. Likewise, the 
two parameters which are relatively constant in the 
phugoid but which vary considerably in the short period, 
angle of attack and pitch rate, had little effect оп 


the phugoid but significantly damped the short period. 
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feel. THE ORBITER 


The Aerodvnamic Design Data Book was obtained from 
mies Lyndon B. Johnson Space Center in order to obtain 
data necessary to calculate the stability derivatives and 
Memacquire the block diagram of the longitudinal control 


system (automatic mode) of the space shuttle orbiter. 


The environment chosen for the study was the area at 
which the Mach number was 1,05. This area was chosen 
because the aerodynamics are undergoing rapid change here, 
EU? would provide a testing challenge for programming 


End solution. 


The stability derivatives were found by plotting 
the force and moment coefficients with respect to velocity 
(U) and angle of attack and then measuring the slope; see 
Figures 8-1 through 8-6. By choosing points close together 
where the slope was rapidly changing, piecewise linearity 
could be assumed between them. The results are shown in 
Table 8-1. The derivatives with respect to de were also 
taken from the publication; the graph of ACp, -vs- dy 
is shown, as an example, in Figure 8-7. The 6 derivatives 


menli sted in Table 8-2. 
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0.0500 
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In cases where U was the independent variable, angle 
of attack was assumed to be 10%; where angle of attack 
was the independent variable, U was 1016.9 fps (Mach = 
1.05). The de derivatives were based upon an angle of 


баск оғ 109, 


The point in the design reentry trajectory where the 
Mach number is 1.05 occurs at approximately 62000 feet in 
altitude. Also, at this point, the load factor n, = 


1.09, and the air density Р = 2.0325E-4 slues/ft?. 


Other constants which are required for computations 


асе: 


wing area = $ = 2690 ft? 
wing DPI с = 39.56 Er 
moment of inertia = Ly о slug-ft^ 


gravity = с = 32.073 ft/sec? (corrected 


tor Altitude and latitude) 


Е 2> @ 9.0 deg 


mass = m 5839 slugs 


For the condition of Mach = 1.05 and angle of attack 
= 10°, the stability derivatives take on the following 
values: 
eV? SE C 

m (ztm- >) = - 004410 


> F: 
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ES Е-5 
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Entering these values into the plant matrix (see page 


54) yields 
-0,04410 0,02849 0.0 -32.073 


| ) -0.09845 -0.01635 101639 0:0 
ӨШ = 


-4 ,683E-4 -9.730E-6 -0.09877 000 


0.0 0.0 1.0 0.0 


-0. 144219 
| -0.78489 
{в} = 
-0,01476 
0.0 ) 


Entering BASMAT with [4] yields the results shown in 
Figure 8-8. Using the GRAPH program developed in Chapter 
VI (see page 71), the time response to an initial condi- 
ШИ 07 мо = 10 fps was plotted. The results are shown 
Meni cure 8-9, The vehicle is clearly unstable in the 


open-loop configuration. 


A look at the eigenvalues printed out in Figure 8-8 
showed that the vhugoid had disappeared, being replaced 
by two real roots, one of which was positive. Although 
the negative (stable) root was larger, the instability 


became dominant very quickly. 


The eigenvalues were not widely separated, having 
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periods of 75.1, 61.2, and 39.7 seconds, respectively, 
so the short period was effectively missing as well. 
EE uouzn its theoretical damping ratio was 0.4337, it 
was completely supressed by the instability of the 


positive real root. 


It was odvious that some form of feedback network 
would be necessary so that the stability of the vehicle 
ТІП be more in line with that of conventional aircraft. 
And, since the orbiter was designed to fly at the very 
extremes of atmosoheric conditions, the system would 
necessarily depend upon many varameters, among them 
dynamic vressure, Mach number, true airsveed, load factor 
Euch is a function of C), pitch angle and rate, angle 
of attack, and elevon deflection angle. Other factors 
could have been used, e.2., vertical velocity, descent 
angle, to name a few, but the ones chosen appeared to be 
those most easily obtainable and from which other 
necessary information could be computed. The end result 


of the engineering design is shown in Figure 8-10. 


In order КО СУТ үе tne equations of motion as they 
were modified by the various gains, feedbacks, lead-lag 
elements, and limits, the Continuous System Modeling 
Program (version III), as incorporated on the IBM360, 


Mas utilized. 
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The system variables were assigned names as follows: 
N, command = GCMD 
true airspeed = Ul 
N, “E = NZ 


Я = ТНЕТА1 


gorizontal velocity 
о = X1 


vertical velocity 
perturbation = X2 


КАС rate perturbation = X3 
pitch angle perturbation = X4 
elevon deflection 
pescurbacaon = X8 
The invut signal to the first lead-lag element was 
called X5IN and its output called X5. The input to the 


Second lead-lag element was called XÓIN and its output X6. 


The signal entering the % LIMIT was called X8PRE1 
and the limited signal X8PRE2; the signal entering the 
be INIT was called X8PRE3 and its output called X8DOT; 
the integration of X8DOT was X8. The signal X8 was fed 


back through a las element whose output was labeled X7. 


The CSNP IIl program is divided into three segments, 
ШҮШТАТ,, DYNAMIC, and TERMINAL. The INITIAL part is for 


defining parameters whose value will not change for the 
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duration of the program. The DYNAMIC segment contains 
the equations of motion and other computations. The 
TERMINAL section (optional) contains any final operations 
which are desired. For this analysis, the following 


terms were assumed constant: 


gravity = GRVTY = 32.073 ft/sec? 


Ag = да - 0.0 


а а 

density = REO = 2.0325E-4 sluz/ft2 

Сид + CMADOT = -2.0 
The initial values of two parameters were listed: 

pitch angle = THETA = 0.1571 rad 

I| e psseed - (0 = 1016.9 Tos 
The other constants listed in the program (see Figure 
8-11) are dimensional parameters or mass properties of 
the orbiter. AREA refers to wetted wing area; ISUBY 
(Ty) is the mass moment of inertia about the pitch axis; 


CHORD is the reference wing chord length; MASS is the 


mass of the vehicle with the payload out. 


The FUNCTION statements enable the program to find 
a dependent variable knowing the value of the independent 
variable, much like reading a graph. The interpolation 
between the points listed in the statement may be linear 


or quadratic. Should an independent variable value fall 
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outside the range of points listed, the program will 


extrapolate a value based on the slope of the function 


at the nearest point listed. 


ЕО ТОМ statements are entered with an AFGEN 


command (for linear interpolation). Since linear inter- 


polation was used in this 


analysis, it was necessary to 


insure that the dependent variable varied linearly 


between the points chosen 


#1; CMDE 
E2: ODDE 
po CEDE 
F4: GCMD 
Pose CUA 
org 
NIU 
pe 
F9: CDA 
Eon CD) 
СИУ МО 
Е12: СМ = 
Е13: СМА 
F14: CMQ 


The functions were 
= Cus 
Со 
= (Di 
= 4 Command 


A 


См. 
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These coefficients were in turn used to compute the 
stability derivatives. As the values of the independent 
variables changed, so did the dependent variables and so 
aso did the stability derivatives. Іп this manner, the 
stability derivatives were modified due to changing 
merent conditions, and thus an important facet of the 
model was maintained. The stability derivatives were 
then used to write the equations of motion (the X1DOT 


through XY4DOT equations). 


The modeling of the system was completed by writing 
four other equations as the biock diagram signified. 
Х5 апа Х6 were written simply as lead-lag (LEDLAG) 
mame tions of X5IN and X6IN, respectively, where the 
parameters listed were the coefficients of s in the 


mamerator and denominator. 


The term X7 was found by writing the lag compensation 
com which it came as a differential equation. Since 
output divided by invut equals the transfer function, 

X7 Бао 
A SS 
ОСИБ 0) = 1,5(X8) 


Taking the lnverse Laplace transform, 
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X7 Eu o cM. 5 CX 


and X7 = / X? 


Writing this in the CSMP III language, 
ЕСО. 0, ХОТ) 
Pere the 90.0 refered to the initial condition. 


The equation to find X8 was somewhat more complex. 
The sisnal was first the sum of GDQ*X6 + Х? - XBPRE1; 
this was then limited to -35/+209, comoared to X8, and 
poesdiffterence multiplied by 20 (here, in radians). 
After limiting this signal to -20/+200 per second, it was 


iicegrated to find X8. 


The eight equations were solved simultaneously using 
the STIFF method of integration, and, after converting 
some of the parameters into degrees, the results were 


printed > 


Шота ОР seen from the block diagram that the forcing 
ШИС стол 15 the command load factor, CCND. This value, 
plotted as FUNCTION F4, was available from the Aero- 


dynamic Design Data Book and was programmed accordingly. 
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The value of CMADOT was not available from the 
Above publication. To arrive at the value of -2.0, 
several aircraft were studied, considering wing plan- 
ШӘП, аіусгатт size, airspeed, and altitude. The larger 
aircraft (C-47 and DC-8) had values of -10 and -6.8, 
respectively, the DC-3 figure being at Mach 0.88 and 
33000 feet. The A-4D, having a delta vlanform, had a 
Cue of -2.1 at 15000 feet and -1.4 ae 35000 feet, 
Beyisidering the size of the orbiter and its altitude 


for this study, the value of -2.0 was deemed appropriate. 


Several equations to compute the load factor on the 
vehicle were considered. The one chosen was derived from 
the lift equation: 


nz = - where L = 0185 


CrdS/ma 


CL*QHAT* AREA/(MASS+CRVTY 


SO, ng 


This equation was chosen as it was most independent of 
the equations of motion and would therefore amplify to 
a less extent any errors due to integration round-off or 
to other sources during the solution of the Seal taneous 


equations. 
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The results are shown in Figure 8-12 (short period) 
and Figure 8-13 (phugoid). Inasmuch as both modes were 
evidently stable, the steady state magnitudes were 


extremely high. Several causes must be considered. 


One, the equations of motion were written as per- 
turbed motions. This was necessary in order to justify 
their linearization. The parameter values were treated 
as perturbations throughout the prozram, being added to 
the large scale motion where necessary for computation 
en the case ot THETA! and U1) and for output plots. 
Mets initial conditions were collectively zero. The 
вета! conditions for the other parameters, X5 through 
X7, were not zero, however, and a method for incorpora- 
ting (or even deriving) their initial value was not 
available. The initial condition for X8 was available 
from the Data Book and was so used. t 1s not known 
the Jack of initial conditions for these parameters 
caused the large excursions in the response; it is 
doubtful that this was the case because the steady state 


values should still have been reasonable; they were not. 


Two, the simulation was begun at some voint other 


than t=0 in the trajectory. That is, the vehicle had 
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already entered the atmosphere, was flying aerodynamically, 
being commanded, responding, etc., in short, was in a 
dynamic environment when this program started. Unless 

ШІБ parameters existing at t, for this simulation can be 
precisely ascertained, it is doubtful that meaningful 


results can be expected. 


Ше ШПІ5 End, angle of attack Corresponding to the 


commanded load factor at t. was calculated: 


М] 


г 
a 


= 1,09 = CL(3pV2)S/me 


КО О TOF CU, yielded 0.722. Since the angle of attack 
iO "s" was 109 with a Cy; of 0.544, the angle of 
attack at 1.09 "с" was found to be ОО ИР е Т) 


8 13,279 (0,2316 гай). 


КОШ ТТЕШ рісі гасе would affect the X2 coefficient 
ши Не Х РОТ equation and the Xi coefficient in the X2DOT 
equation, so it was computed using n, - QoUo/g + 1.0 


and found to be 0.0344 rad/seco 


Poth of these modifications were made to the program 


with negligible effect on the results. 


To analyze the results as they stand: 


The angle of attack first exceeded the limits of the 
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FUNCTION statements (requiring interpolation of the force 
coefficients well outside the linear range of the functions) 
at approximately 2 seconds, oscillated in and out of the 
range, and left the range for good at approximately 10 
seconds. All stability derivatives dependent upon angle 

of attack must be considered suspect after this time. 

The angle of attack was computed by dividing the vertical 
velocity perturbation by the horizontal velocity and adding 


із to the initial angle of attack. 


The vertical velocity response (positive is downward) 
reached an extremely nigh value very quickly, and if applied 
to altitude, would have the orbiter impacting the ground 
before the simulation period was half over. Such a rate 
of descent would also cause a rapid change in air density, 
which was orizinally assumed constant, leading to very 
 Гегепо values of those stability derivatives proportional 
THO. The assumption of a constant RHO should be valid 


mem сле three minute-duration of the simulation. 


The pitch rate magnitudes are not unusual, perhaps a 
bit high, but the steady state pitch angle is much too 
ШЕ!) The steady state values of 58° pitch angle and 71° 
angle of attack indicate a glide angle of 13°, also much 


poo high. 


188 





The elevon deflection, going to the trailing edge up 
limit in just over ten seconds, indicated that the vehicle 
sensed an insufficient nose up vitching moment with a 
GCMD of 1.09 and was trying vainly to create that moment. 
The result of the full trailing up deflection was instead 


a stalling angle of attack and a high rate of descent. 


The actual load factor experienced by the vehicle in 
the simulation was only slightly higher than expected; the 
large elevon deflection was undoubtedly the cause of the 
initial ravid rise in this parameter., However, since the 
BeoOgram computed nz as a function of p , the steady state 


value was meaningless. 


Taking the above observations into consideration, the 
program appeared to be consistent. Many possible reasons 
for the magnitude discrevancies were examined, but no 
changes in the vrogram made more than a modicum of change 
ln Lhe results. It appeared that the error was basic in 
Nature and not related to the model as presented or to 
the CSMP program. The most likely cause is, as was 
mentioned earlier, the fact that tne simulation is not 


Semmeicing at t=0. 
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